!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
recursive subroutine hilbert_space_filling_2d (base, current, start_pos, end_pos, start_point, order, var_regrid_flag)
!==================================================================================================================================|

   use parameters
   use block_define
   implicit none

   integer(4), intent(in) :: order, var_regrid_flag
   integer(4), dimension(2), intent(in) :: start_pos, end_pos, start_point
   type(block), pointer, intent(in) :: base
   type(block), pointer, intent(inout) :: current

   integer(4) :: order_new
   integer(4), dimension(2) :: end_pos_r, start_pos_r, start_point_new

!----------------------------------------------------------------------------------------------------------------------------------|
   end_pos_r = 1 - end_pos
   start_pos_r = 1 - start_pos

   if (order .eq. 1) then
      start_point_new = start_point

      if (start_point_new(1) .le. base%nc_x .and. start_point_new(2) .le. base%nc_y) then
         if (var_regrid_flag .eq. cpu_id) then
            current%next => base%child(start_point(1), start_point(2))%pblock
            current%next%last => current
         endif
         current%framework_next => base%child(start_point(1), start_point(2))%pblock
         current%framework_next%framework_last => current
         current => current%framework_next
         current%start_pos = start_pos
         current%end_pos = end_pos_r
         if (start_pos(1) .eq. 0) then
            current%start_point(1) = 1
         elseif (start_pos(1) .eq. 1) then
            current%start_point(1) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
         if (start_pos(2) .eq. 0) then
            current%start_point(2) = 1
         elseif (start_pos(2) .eq. 1) then
            current%start_point(2) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
      endif

      start_point_new = start_point_new + end_pos_r - start_pos
      if (start_point_new(1) .le. base%nc_x .and. start_point_new(2) .le. base%nc_y) then
         if (var_regrid_flag .eq. cpu_id) then
            current%next => base%child(start_point_new(1), start_point_new(2))%pblock
            current%next%last => current
         endif
         current%framework_next => base%child(start_point_new(1), start_point_new(2))%pblock
         current%framework_next%framework_last => current
         current => current%framework_next
         current%start_pos = start_pos
         current%end_pos = end_pos
         if (start_pos(1) .eq. 0) then
            current%start_point(1) = 1
         elseif (start_pos(1) .eq. 1) then
            current%start_point(1) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
         if (start_pos(2) .eq. 0) then
            current%start_point(2) = 1
         elseif (start_pos(2) .eq. 1) then
            current%start_point(2) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
      endif

      start_point_new = start_point_new + end_pos - start_pos
      if (start_point_new(1) .le. base%nc_x .and. start_point_new(2) .le. base%nc_y) then
         if (var_regrid_flag .eq. cpu_id) then
            current%next => base%child(start_point_new(1), start_point_new(2))%pblock
            current%next%last => current
         endif
         current%framework_next => base%child(start_point_new(1), start_point_new(2))%pblock
         current%framework_next%framework_last => current
         current => current%framework_next
         current%start_pos = start_pos
         current%end_pos = end_pos
         if (start_pos(1) .eq. 0) then
            current%start_point(1) = 1
         elseif (start_pos(1) .eq. 1) then
            current%start_point(1) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
         if (start_pos(2) .eq. 0) then
            current%start_point(2) = 1
         elseif (start_pos(2) .eq. 1) then
            current%start_point(2) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
      endif

      start_point_new = start_point_new + end_pos - start_pos_r
      if (start_point_new(1) .le. base%nc_x .and. start_point_new(2) .le. base%nc_y) then
         if (var_regrid_flag .eq. cpu_id) then
            current%next => base%child(start_point_new(1), start_point_new(2))%pblock
            current%next%last => current
         endif
         current%framework_next => base%child(start_point_new(1), start_point_new(2))%pblock
         current%framework_next%framework_last => current
         current => current%framework_next
         current%start_pos = start_pos_r
         current%end_pos = end_pos
         if (start_pos_r(1) .eq. 0) then
            current%start_point(1) = 1
         elseif (start_pos_r(1) .eq. 1) then
            current%start_point(1) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
         if (start_pos_r(2) .eq. 0) then
            current%start_point(2) = 1
         elseif (start_pos_r(2) .eq. 1) then
            current%start_point(2) = ratio
         else
            write (*, *) 'Hilbert curve error from subroutine hilbert_space_filling_2d'
            stop
         endif
      endif
   else
      order_new = order - 1

      start_point_new = start_point
      call hilbert_space_filling_2d (base, current, start_pos, end_pos_r, start_point_new, order_new, var_regrid_flag)

      start_point_new(1) = start_point_new(1) + (end_pos_r(1) - start_pos(1)) * 2 ** order_new
      start_point_new(2) = start_point_new(2) + (end_pos_r(2) - start_pos(2)) * 2 ** order_new
      call hilbert_space_filling_2d (base, current, start_pos, end_pos, start_point_new, order_new, var_regrid_flag)

      start_point_new(1) = start_point_new(1) + (end_pos(1) - start_pos(1)) * 2 ** order_new
      start_point_new(2) = start_point_new(2) + (end_pos(2) - start_pos(2)) * 2 ** order_new
      call hilbert_space_filling_2d (base, current, start_pos, end_pos, start_point_new, order_new, var_regrid_flag)

      start_point_new(1) = start_point_new(1) + (end_pos(1) - start_pos(1)) * (2 ** order_new - 1)
      start_point_new(2) = start_point_new(2) + (end_pos(2) - start_pos(2)) * (2 ** order_new - 1)
      start_point_new(1) = start_point_new(1) + (end_pos(1) - start_pos_r(1))
      start_point_new(2) = start_point_new(2) + (end_pos(2) - start_pos_r(2))
      call hilbert_space_filling_2d (base, current, start_pos_r, end_pos, start_point_new, order_new, var_regrid_flag)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine hilbert_space_filling_2d
